Multiple dark-bright solitons in atomic Bose-Einstein condensates 
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We present experimental results and a systematic theoretical analysis of dark-bright soliton interactions and 
multiple-dark-bright soliton complexes in atomic two-component Bose-Einstein condensates. We study ana- 
lytically the interactions between two-dark-bright solitons in a homogeneous condensate and, then, extend our 
considerations to the presence of the trap. An effective equation of motion is derived for the dark-bright soliton 
center and the existence and stability of stationaiy two-dark-bright soliton states is illustrated (with the bright 
components being either in- or out-of-phase). The equation of motion provides the characteristic oscillation 
frequencies of the solitons, in good agreement with the eigenfrequencies of the anomalous modes of the system. 
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I. INTRODUCTION 

Over the past few years, the macroscopic nonlinear struc- 
tures that can be supported in atomic Bose-Einstein conden- 
sates (BECs) have been a topic of intense investigation (see, 
e.g., Refs. lol-H] for reviews in this topic). The first experi- 
mental efforts to identify the predominant nonlinear structure 
in BECs with repulsive interatomic interactions, namely the 
dark soliton, were initiated over a decade ago IH-lst. However, 
these efforts suffered from a number of instabilities arising 
due to dimensionality and/or temperature effects. More re- 
cently, a new generation of relevant experiments has emerged, 
that has enabled the overcoming (or quantification) of some 
of the above limitations. The latter works have finally enabled 
the realization of oscillating, and even interacting, robust dark 
solitons in atomic BECs. This has been achieved by means of 
various techniques, including phase-imprinting/density engi- 
neering 1 10- 12], matter-wave interferen ce |13l fl4ll . or drag- 
ging localized defects through the BECs II15I1 . 

Atomic dark solitons may also exist in multi-component 
condensates, where they are coupled with other nonlin- 
ear macroscopic structures lUl HI Si- Of particular inter- 
est are dark-bright (DB) solitons that are supported in two- 
component 1 16] and spinor (vf\ condensates. Such struc- 
tures, are frequently called "symbiotic" solitons, as the bright- 
soliton component (which is generically supported in BECs 
with attractive interactions fs*]) may only exist due to the inter- 
species interaction with the dark-soliton component. Dark- 
bright solitons have also attracted much attention in other 
contexts, such as nonlinear optics |18] and mathematical 
physics ifTgll . In fact, DB-soliton states were first observed 
in optics experiments, where they were created in photore- 
fractive crystals |.2ftl. while their interactions were partially 
monitored in Ref. ll2lll . In the physics of BECs, robust DB- 
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solitons were first observed in the experiment of Ref. [To'] by 
means of a phase-imprinting method, and more recently in 
Refs. 1I22I - I24I1 by means of the counterflow of the two BEC 
components. The above efforts led to a renewed interest in 
theoretical aspects of this theme: this way, DB-soliton inter- 
actions were studied from the viewpoint of the integrable sys- 
tems theory in Ref. Izsl], DB-soliton dynamics were investi- 
gated numerically in Ref. [26], while DB-solitons in discrete 
settings were recently analyzed in Ref. |27]. Furthermore, 
higher-dimensional generalizations — namely, vortex-bright- 
soliton structures — were recently studied as well [28]. 

Our aim in the present work is to study multiple-DB soli- 
tons in two-component BECs confined in harmonic traps. 
First, we present our experimental results, based on the coun- 
terflow of two rubidium condensate species [22-24], which 
demonstrate the existence (and indicate the robustness) of 
such multiple-DB-soliton states. Motivated by the experimen- 
tal observations, we then proceed to analyze the interactions 
of DB solitons, first in the case of a homogeneous system and, 
next, in the presence of the trap. Our analytical approximation 
relies on a Hamiltonian perturbation theory, which leads to an 
equation of motion of the centers of DB-soliton interacting 
pairs. Employing this equation of motion, we demonstrate the 
existence of stationary two- and three-DB-soliton states, find 
semi-analytically the equilibrium distance of the constituent 
solitons, as well as the oscillation frequencies around these 
equilibria. The oscillation frequencies correspond to the char- 
acteristic anomalous modes' eigenfrequencies that we com- 
pute via a Bogoliubov-de Gennes (BdG) analysis. This way, 
we are able to quantify the properties of stationary multiple- 
DB-solitons in harmonically confined two-component BECs, 
and provide analytical results for their in- and out-of-phase 
motions. 

The paper is organized as follows. In Section II, we provide 
our experimental results. In Section III we describe our the- 
oretical setup and present the DB-soliton states. Section IV 
is devoted to the study of the interactions of two DB-solitons, 
while Section V contains the results for multiple DB-solitons 
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FIG. 1: (Color online) Experimental images indicating DB-soliton 
clustering in a two-component BEC. The upper cloud in each image 
(and red curve in inset) shows atoms in the |2, — 2) state, while the 
lower cloud (black curve) shows atoms in the |1, — 1) state. Prior to 
imaging, the two components are overlapped in trap for 5 sec. Insets 
show integrated cross sections of the boxed region. For details see 
text. 



in the trap. Finally, in Section VI we summarize our findings 
and discuss future challenges. 



magnetic gradient of 10.4 mG/cm is applied along the elon- 
gated axis of the BEC, inducing counterflow of the two com- 
ponents. As a result, a dense modulation instability pattern 
arises. Once the pattern has fully developed, the gradient is 
turned off. During the subsequent in-trap evolution, the ini- 
tially very regular pattern becomes irregular, and over a time 
scale of several seconds displays the dynamics of interacting 
solitons originating from the modulational instability. Images 
taken in several experimental runs 5 sec after the switch-off of 
the gradient are presented in Fig. [T] in which the upper cloud 
(and the red curves in the insets) in each image shows atoms 
in the |2,— 2) state after 7 ms of free expansion, while the 
lower cloud (and black curves in the insets) shows the atoms 
in the |1, — 1) state after 8 ms of expansion. This difference in 
expansion time just serves to separate the two states vertically 
for imaging. 

An intriguing observation is the frequent formation of large 
gaps in one component (which constitutes the component sup- 
porting the dark-solitons) that are filled by bright-solitons in 
the other component. Interestingly, these gaps are structured 
by small, periodic density bumps, indicating that these regions 
are composed of merged solitons. Some of these features are 
marked by the boxed regions in Fig. [T] with corresponding 
cross sections shown as insets. We clearly observe clusters of 
two- and three-merged solitons [see Fig.[Tla-c)], and also have 
some indications of clusters composed of four- to five-solitons 
— see Fig.fTJd, e). 

While our destructive imaging technique does not allow us 
to analyze the dynamics and lifetime of the clusters in detail, 
the occurrence of large DB-soliton clusters strongly supports 
the theoretical part of our work that we will present below: in 
fact, we will study analytically the interaction between two- 
DB solitons, and we will demonstrate the existence of two- 
and multiple-DB stationary states, resembling the ones ob- 
served in the experiment. Furthermore, we will study the sta- 
bility of these states and discuss their dynamics in the pres- 
ence of the harmonic trap. 



II. EXPERIMENTAL RESULTS - MOTIVATION 

Since our scope is the study of multiple-DB-solitons in 
atomic BECs, we start by presenting some experimental re- 
sults, which showcase the existence of such structures. These 
results, apart from being interesting in their own right — as 
they demonstrate the formation of DB-soliton clusters in two- 
component BECs — provide the motivation for a systematic 
analysis of muItiple-DB-solitons, which will be presented in 
the following sections. 

Our soliton generation scheme is based on a counterflow 
induced modulational instability details of which have been 
described in Refs. [22, 23]. Briefly, we start with a BEC of 
about 800, 000 ^^Rb atoms in the |F, mp) = |1, -1) hyper- 
fine state. The atoms are confined in an elongated optical 
dipole trap with measured trap frequencies of 27rx{1.5, 140, 
178} Hz. About half the atoms are then transferred to the 
1 2,— 2) state with a brief microwave sweep, thus producing 
a weakly miscible two-component mixture. Subsequently, a 



III. MODEL AND THEORETICAL SETUP 

A. Coupled GPEs and dark-bright solitons 

Following the experimental observations of the previous 
section, we consider a two-component elongated (along the 
x-direction) BEC, composed of two different hyperfine states 
of rubidium. As is the case of the experiment, we consider 
a highly anisotropic trap, with the longitudinal and transverse 
trapping frequencies such that uj^ ^ In the framework 
of the mean-field theory, the dynamics of this two-component 
BEC can be described by the following system of two coupled 
GPEs liii: 

ihdtipj=\-—dl'il^j + V{x) - ^lJ +^5jfc|V'fcNV'j- (1) 

Here, ipj{x,t) (j = 1,2) denote the mean-field wave func- 
tions of the two components (normalized to the numbers of 



3 



atoms Nj = \%l)j\'^dx), m is the atomic mass, /ij are 

the chemical potentials, and V{x) represents the external har- 
monic trapping potential, V{x) ~ (l/2)mf2^x^ where Vl = 
lUx/ijo^- In addition, gjk = 2huj±ajk are the effective ID cou- 
pling constants, ajk denote the three s-wave scattering lengths 
(note that ai2 = 021) accounting for collisions between atoms 
belonging to the same (ajj) or different iajk,j 7^ k) species. 
In the case of the hyperfine states |1, —1) and |2, —2) of ^^Rb 
considered in the previous section, the scattering lengths take 
the values an — 100. 4ao, 012 = 98.98ao and 022 = 98.98ao 
(where oq is the Bohr radius) ||22ll23ll . Thus, we may safely 
use the approximation that all scattering lengths take the same 
value, say w a [38]. To this end, measuring the densities 
length, time and energy in units of 2a, a± — ^ h/tuj^, 
wj^ and fio;^, respectively, we may reduce the system of 
Eqs. ^ into the following dimensionless form. 



+ (IV'jf + |^3-,f -MjOV-,, J = 1,2. (2) 

Below, we will consider a situation where the component 
characterized by the wavefunction (1/12) supports a single- 
or a multiple-dark (bright) soliton state, and the respective 
chemical potentials will be such that /ii > /i2. Note that 
the external potential in Eqs. (|2]l takes the form V{x) = 
{1 / 2)0,'^ x"^ , where fl = uJx/oj^ ^ 1 is the normalized trap 
strength. 

We assume that a single- or a multiple-dark-soliton state is 
on top of a Thomas-Fermi (TF) cloud with density It/itfP = 
III — V{x); this way, the density in Eqs. (|2]l is substi- 
tuted as iV'iP iV'TFplV'iP- Furthermore, we introduce the 
transformations t — >■ fiit, x — y/JIix, \tp2\'^ /^r^l^2p, 
and cast Eqs. (|2]i into the following form: 



l)V'i=i?i, (3) 



(|V'1|' + IV'2P-M)V'2 =i?2, (4) 



where jl = fi2/l^i, while 

i?i = (2/i?)-^ [2(l-|V^i|2)y(a;)Vi 

R2 = 11^^(1 -\iji\^)V{x)ij2], 



V'{x)d,^i 



(5) 



with V'{x) = dV/dx. Equations ©-dUl can be viewed 
as a system of two coupled perturbed nonlinear Schrodinger 
(NLS) equations, with perturbations given by Eqs. (|5}. In the 
absence of the trap (i.e., for O = 0), the perturbations vanish 
and Eqs. ©-(Illi actually constitute the completely integrable 
Manakov system ll29ll . This system conserves, among other 
quantities, the Hamiltonian (total energy). 
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as well as the total number of atoms, N = Ni + N2 = 
Y^^=i f-^ \ipj\'^dx; additionally, the number of atoms of 
each component, iVi and A'2, is separately conserved. 

Considering the boundary conditions IV'iP 1 and 

|V'2p -J> as |a;| 00, the NLS Eqs. ©-Q possess an 
exact analytical single-DB soliton solution of the following 
form (see, e.g., Ref. 116(1 ): 

ipi{x,t) = cos(j)ta,nh[D{x — xo{t)] + isincj), (7) 
ip2{x,t) = r]sech[D{x — xo{t)]exp[ikx + i6{t)] , (8) 

where is the dark soliton's phase angle, cos0 and 77 rep- 
resent the amplitudes of the dark and bright solitons, D and 
xo{t) denote the width and the center of the DB soliton, while 
k = Dtancf) ~ const, and 9(t) are the wavenumber and 
phase of the bright soliton, respectively. The above parame- 
ters of the single DB-soliton are connected through the fol- 
lowing equations: 



Xq 

0{t) 



cos 



D tan ( 
1 



-AD' -k')t+{fL-l)t, 



(9) 
(10) 
(11) 



where io = dxQ jdt is the DB soliton velocity. Below, we will 
mainly focus on stationary solutions, characterized by a dark 
soliton's phase angle = [in this case, the bright soliton 
component is stationary as well — see Eq. (fTOliI: nevertheless, 
we will also consider the near-equilibrium motion of DB soli- 
tons, characterized by (/) « 0. 

To describe a two-DB-soliton state (for £7 = 0) composed 
by a pair of two equal-amplitude single DB solitons traveling 
in opposite directions, we will use the following ansatz: 

ipilxjt) — (cos (/) tanh X_ + i sin 0) 

X (cos (/i tanh — i sin 0) , (12) 

M^A) = ?7SGchX_e^[+'=-+''W+('^-i)*l 

+ rysechX+e*[-^-+«W+('^-l)*le^^^ (13) 

where X± — D {x zt Xo{t)), 2xq is the relative distance be- 
tween the two solitons, and A9 is the relative phase between 
the two bright solitons (assumed to be constant); below we 
will consider both the out-of-phase case, with A9 = tt, as 
well as the in-phase case, corresponding to A6 ~ 0. 

Notice that in either cases of single- or multiple-DB- 
solitons, the number of atoms of the bright soliton, N2, may 
be used to connect the amplitude 77 of the bright soliton(s), 
the chemical potential /ii of the dark soliton(s) component, as 
well as the width D of the DB soliton. In particular, in the 
case of a single-DB-soliton, one finds that N2 — 2rf',JJIi/ D 
[for the variables appearing in Eqs. (|2)], while for the case of 
a two-DB-soliton state (with well-separated solitons) the rele- 
vant result is approximately twice as large, namely: 



(6) 



N2 



D 



(14) 
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B. Stationary states and their excitation spectrum 

In our numerical computations, we will initially obtain 
(by means of a fixed-point algorithm) stationary solutions 
of Eqs. (|2|i, in the form ipi{x,t) = u{x) and tp2{x,t) — 
v{x), and then consider their linear stability. This is numer- 
ically studied via the well-known BdG analysis (see, e.g., 
Refs. iQl, 12, 01), upon introducing the following ansatz into 
Eqs. Q, 



ipi{x,t) = u{x)+e a{x)e^* + b* {x)e 



ip2{x,t) = v{x)+e c{x)e^^ + d* {x)e 



(15) 
(16) 



The resulting equations are linearized (keeping only terms 
of order of the small parameter e), and the ensuing eigen- 
value problem for eigenmodes {a{x),b{x),c{x),d{x)} and 
eigenvalues A = Ar + iXi is solved [note that the asterisk 
in Eqs. ([TSll-dTSIl denotes complex conjugation]. In the case 
of a single DB soliton, the excitation spectrum can be well- 
understood in both cases, corresponding to the absence and 
the presence of the harmonic trap, using the following argu- 
ments. 

First, in the absence of the trap, the system of Eqs. ^ fea- 
tures not only a U{1) (phase) invariance in each of the com- 
ponents but also a translational invariance; thus, the system 
has three pairs of eigenvalues (each associated with one of the 
above symmetries) at the origin of the spectral plane (Ar, Ai). 
In this case, the phonon band (associated with the continuous 
spectrum of the problem) covers the entire imaginary axis of 
the spectral plane. 

Second, in the presence of the trap, the single DB soli- 
ton "lives" on the background of the confined ground state, 
i.e., {tpi'jl^2} = {'0TF,O} (as discussed above). It is well- 
known fv, that the harmonic potential introduces a dis- 
crete (point) BdG spectrum for this spatially confined ground 
state. In addition to that, the translational invariance of the 
unconfined system is broken and, due to the presence of the 
DB soliton, a single eigenvalue A'^^'^^ emerges. The respec- 
tive (negative energy) eigenmode is the so-called anomalous 
mode (AM), while the associated eigenvalue A*^^'^-' is directly 
connected with the oscillation frequency of the DB soliton in 
the harmonic trap, similarly to the case of a dark soliton in 
one-component BECs |30]. In fact, the imaginary part of the 

eigenvalue A*^^'^^ reads A^^'^'' = Wosc, where Wqsc is the os- 
cillation frequency of the single DB soliton, given by 11611 : 



1 



X 



Xo 
Xo 
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/^J.l 



1 
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(17) 
(18) 



The above results are illustrated in Fig. |2] where a typi- 
cal example of a stationary single DB-soliton state is depicted 
(top panel); additionally, the eigenvalues A^ characterizing the 
excitation (BdG) spectra of such stationary states, are shown 
as functions of the chemical potentials fii and fj.2 in the mid- 
dle and bottom panels of the figure, respectively. As observed 




FIG. 2: (Color online) The top panel depicts the stationary solution 
for a single DB-soliton for fii = 3/2, ~ 1, and Q, = 0.1. The 
bright (dark) components are shown by the dashed green (solid blue) 
lines. The middle (bottom) panel shows the normalized imaginary 
part Ai/f2 of the eigenvalues for the single DB-soliton as a function 
of fii (jj,2) for jj,2 ~ 1 = 3/2). The (red) dashed line, depicts 
the analytical prediction of Ref. fl?] for the DB-soliton oscillation 
frequency [cf. Eq. ( I17H . providing an excellent approximation to the 
anomalous mode eigenfrequency. 



in these two bottom panels, there exist two types of spectral 
lines, namely "slowly-varying" ones (analogous to ones that 
are present in the spectrum of a dark soliton in one-component 
BECs fisll ) and "fast-varying" ones due to the presence of the 
second (bright-soliton) component. The latter, as was pointed 
out also in the recent work of Ref. ll24ll may, in fact, col- 
lide with the internal anomalous mode of the DB soliton and 
give rise to instability quartets which are barely discernible in 
Fig. |2](see, e.g., the bottom panel for /i2 > 1.4 where a merger 
of eigenvalues occurs). Generally, however, it is found that the 
analytical prediction (red dashed line) is excellent in capturing 
the relevant anomalous mode pertaining to the DB-soliton os- 
cillation. 

The above discussion sets the stage for the presentation of 
our results for multiple DB-soliton states. 
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IV. INTERACTION BETWEEN TWO DARK-BRIGHT 
SOLITONS 

We start by considering the case where the external trap is 
absent, i.e., for = 0. To study the interaction of two iden- 
tical DB solitons, as described in the ansatz of Eqs. (fTSji-lfTsTl, 
we will employ the Hamiltonian approach in the framework 
of the adiabatic approximation of the perturbation theory for 
matter- wave solitons (see, e.g., Refs. In particular, we 

assume that the approximate two-DB-soliton state features an 
adiabatic evolution due to a weak mutual interaction between 
the constituent solitons and, thus, the DB soliton parameters 
become slowly-varying unknown functions of time t. Thus, 
(j) — ^ D -> D{t) and, as a result, Eqs. (|9]l-(fT0li become: 



xo{t) = D{t) tan (j){t), 



(19) 
(20) 

where we have used Eq. ( fT4] l. The evolution of the parameters 
(f>{t), D{t) and xo{t) can then be found by means of the evo- 
lution of the DB soliton energy as follows. First, we substitute 
the ansatz (fT2]i-(fT3]l into Eq. (|6]l and perform the integrations 
under the assumption that the soliton velocity is sufficiently 
small, such that cos(fcx) « 1 (and sin(fca;) ~ 0). Then, we 
further simplify the result assuming that the solitons are well- 
separated, i.e., their relative distance is xq ^ 1. This way, we 
find that the total energy of the system assumes the following 
form: 

E = 2Ei + i?DD + Ebb + 2£'db, (21) 
where Ei is the energy of a single DB soliton, namely. 



El 



2 (m - 1) 



D 



D 



(22) 



while the remaining terms account for the interaction between 
the two DB solitons. In particular, -Edd, ^'BB, and E^b de- 
note, respectively, the interaction energy between the two dark 
solitons, the two bright ones, and the interaction energy be- 
tween the dark soliton of one component and the bright one 
in the other component. The above interaction energies are 
given by the following (approximate) expressions: 



E, 



DD 



16 COS^ I 



D cos^ (j> + D + 2(cos2 (j) - D'^)xq 



3 + 4 cos^ ( 
3D 



(23) 



2D {D (1 - Dxo) - k^xo) + Dx 



X 



XD {2Dxo - 1) (1 + 2 cos^ A9) 



-iDxo 



(24) 



E^b = -4xcos2(/)cosA6'e-2^"=« 



xcos 



16 



cos^ - 16132:0 + 8 



e-*^^", (25) 



where terms of order 0{e~^^^°) and higher have been ne- 
glected (nevertheless, it has been checked that their contri- 
bution does not alter the main results that will be presented 
below). 

Having determined the two-DB-soliton energy [up to order 
0(e^^^^")], we can find the evolution of the soliton param- 
eters from the energy conservation, dE/dt = 0. In fact, we 
focus on the case of low-velocity, almost black solitons (with 
D{t) « and cos(/)(i) w 1), for which energy conservation 
leads to the following nonlinear evolution equation for the DB 
soliton center: 



Xo — -Pint, 
Pint = ^DD + -FbB + 2Fdb • 



(26) 
(27) 



In the above equations, Fi^t is the interaction force between 
the two DB solitons (depending on the soliton coordinate xq), 
which contains the following three distinct contributions: the 
interaction forces Fdd and Fbb between the two dark and two 
bright solitons, respectively, as well as the interaction force 
Fdb of the dark soliton of the one soliton pair with the bright 
sohton of the other pair. These forces have the following form: 



Xo 



1(544- 



352i5^) + 128Do {D^ - l) 



Xo 



X e 



-iDoxo 



(28) 



^^BB = 



Xo 



- 6D, 







iDlxa - 2x 



X DlcosA9e^^^°''° 



x2 

Xo 



^1 + 2cos2 Ad) (-SDoXo + 6) 



771 _ ^ 

-f^DB — 

Xo 

+ ^ 

Xo 



8Do COS AO 



208 

— + 64i:ioa;o 



Doe 



(29) 



(30) 



where D{t) « Do since we are assuming that D{t) « 0. 

The equation of motion for the two-DB-soliton state [cf. 
Eq. ( |26] )1 provides a clear physical picture for the interaction 
between the two DB solitons. In order to better understand 
this result, first we note that the leading order interaction force 
between the bright soliton components is cx exp(— 2_Doa^o) 
and, as a result, it is decaying more slowly for large xq than the 
one between the two dark solitons which is cx exp(— 4Z3oa^o); 
the interaction between dark and bright is also to leading order 
cx exp(— 2£'oa;o). This result is in accordance with earlier 
predictions, where the same dependence of the force over the 
soliton separation was found (see, e.g., Refs. ll3lll and iflllsl 
33] for bright and dark solitons, respectively). 

Let us now consider the role of the bright-soliton compo- 
nent. In its absence, i.e., for x = [cf. Eq. (fT9]l1, it is clear 
that Fbb = ^db = and Eq. ( l26l ) describes the interaction 
between two dark (almost black) solitons; in this case, tak- 
ing into regard that Dq = 1, it can readily be found that the 
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pertinent (repulsive) interaction potential is cx 2exp(— 4a;o), 
which coincides with the result of Ref . ll33ll obtained by means 
of a variational approach (see also a relevant discussion in 
Refs. I4J, Il4tl ). On the other hand, when bright solitons are 
present (i.e., for x 7^ 0), the principal nature of the bright- 
bright-soliton interaction — and also of part of the dark-bright- 
soliton interaction one — depends on the factor cos A6': when 
the relative phase between the bright-soliton components is 
= (A0 = tt), i.e., in the in-phase (out-of-phase) case, 
the interaction is repulsive (attractive). This conclusion stems 
from the fact that the coefficients of the terms oc cos /S.9 are 
positive definite since the parameter Do > 1 for every x > 0. 

According to the above, it is clear that the competition be- 
tween repulsive (for dark solitons) and attractive (for out-of- 
phase bright solitons) forces leads to the emergence of fixed 
points in the equation of motion (|26] | or, in other words, to the 
existence of a stationary two-DB-soliton state ll39tl . Below 
we will demonstrate this effect in detail: we will determine 
the fixed (equilibrium) points, Xcq, as solutions of the tran- 
scendental equation resulting from Eq. (l26t for xq = in 
the out-of-phase case, and subsequently study their stability. 
Nevertheless, before proceeding further, we should mention 
that stationary two-DB-solitons were also found numerically 
and experimentally in Ref. |21] in the context of nonlinear 
optics, but their existence details and stability properties were 
not considered. Additionally, although exact two-DB-soliton 
solutions (as well as A/^-DB-soliton solutions) do exist in the 
Manakov system jzsi [34ll . their compUcated form does not 
allow for a transparent physical description of the relevant dy- 
namics, as provided above. 

Let us now study the stability of the equilibrium points in 
the framework of Eq. (l26T l. Introducing the ansatz Xi^{t) — 
Xeq + 5{t), and linearizing with respect to the small-amplitude 
perturbation 5{t), we derive the following equation: 

5 + uj15^Q, (31) 
where the oscillation frequency loq is given by: 



where the phase-difference between the bright-soliton com- 
ponents is taken to be A6' = tt. Physically speaking, the os- 
cillation frequency represents the internal (out-of-phase) 
motion of the two DB-solitons; in fact, as here we deal with 
the homogeneous case (i.e., in the absence of the trap), the 
in-phase motion of the solitons is associated with the neutral 
translation mode due to the translational invariance of the sys- 
tem (the respective in-phase Goldstone mode has a vanishing 
frequency). 

The above analytical predictions have been compared with 
numerical simulations. First, we have confirmed the existence 
of the stationary two-DB-soliton state (in the out-of-phase 
case); a prototypical example of such a state is shown in the 
top panel of Fig.[3](for fii = 8^2/2 = 3/2). We have also de- 
termined the dependence of the equilibrium soliton positions 
(denoted by in the middle panel of Fig. [3]) and the effective 
frequency [cf. Eq. ( |32] )1 on the chemical potential /i2 of 




^2 



FIG. 3: (Color online) Top panel: A stationary DB-soliton pair: the 
solid (blue) line denotes the two-dark-soliton state [recall that each 
dark soliton is associated with a zero crossing], while the dashed 
(green) line denotes the respective two-bright-soliton state. The 
chemical potentials are /ii = 3/2 and ^2 = 1. Middle panel: the 
equilibrium center of mass xo as a function of the chemical poten- 
tial ^2 (for /^i = 3/2). Stars (in red) denote the analytical pre- 
diction of Eq. ( I26t . while circles (in blue) denote the numerically 
obtained soliton center xq. Bottom panel: the oscillation frequency 
for the out-of-phase motion of the DB-soliton pair as a function of 
/i2 (for /ii = 3/2). Stars (in red) depict the analytical result for ujo 
[cf. Eq. l|32H. while circles (in blue) depict the numerically obtained 
imaginary eigenvalue \i (for the out-of-phase soliton motion) of the 
excitation spectrum. 



the bright soliton component. The respective analytical and 
numerical results are shown in the middle and bottom panels 
of Fig. [3] To obtain the numerical results, we have used a 
(least squares) fitting algorithm to accurately identify the am- 
plitude ?/, inverse width D, and equilibrium center of mass xq 
of the bright component. The numerical findings for xq and 
cjQ (the latter is obtained via a BdG analysis, as the imaginary 
eigenvalue Ai of the stationary two-DB-soliton state) are di- 
rectly compared with the semi-analytical results of Eqs. (l26b 
and ( [32I ), respectively. Taking into account the approximate 
nature of the fitting scheme, we find that there is a very good 
quantitative agreement between the analytical and numerical 
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results (see middle and bottom panels of Fig. O. Notice that 
despite the motion of this eigenvalue through the continuous 
spectrum, no instability is observed in the parametric window 
shown in Fig. [3] 



V. MULTIPLE DARK-BRIGHT SOLITONS IN THE TRAP 

Next, let us consider the case of multiple DB-solitons in the 
presence of the harmonic trap. In the presence of the trap, each 
of the multiple-DB-soliton structures is subject to two forces: 
(a) the restoring force of the trap, Ftr [in the case of a sin- 
gle DB-soliton, this force induces an in-trap oscillation with a 
frequency ujqsc — see Eq. (fTTl)]. and (b) the pairwise interac- 
tion force i^int [cf. Eq. ( |27] )1 with other dark-bright solitons. 
Thus, taking into regard that Ftr = — "^osc^o [16], one may 
write the effective equation of motion for the center a;o of a 
two-DB-soliton state as follows: 



int • 



(33) 



One can thus straightforwardly generalize the above equation 
for TV^-interacting DB-soliton states, similarly to the case of 
multiple dark solitons in one-component BECs ifTsiflA 351. 

It is interesting to observe that, in the presence of the trap, 
the restoring force Ftr can generate equilibrium positions, not 
only for out-of-phase bright solitons (whereby such a state 
could be stationary even without the trap as found above), but 
also for in-phase bright solitons. In the latter case, the re- 
pulsion between both the dark- and the bright-soliton compo- 
nent(s), is balanced by the trap-induced restoring force. In the 
case of two-DB solitons placed at a; = ±xq, the equilibrium 
points, Xcq, can readily be found (as before) as solutions of the 
transcendental equation resulting from Eq. (l3?t for xq = 0, in 
both the in- and out-of-phase cases. To study the stability of 
these equilibrium points in the framework of Eq. (|33) , we may 
again use the ansatz xo{t) — x^q + S{t), and obtain a linear 
equation for the small-amplitude perturbation 5{t), similar to 
that of Eq. ( [3T| i, namely: S + ojfS = 0, where the frequency 
uji is given by. 



(34) 



where luq is given by Eq. ( l32b . Similarly to the case of dark 
solitons in one-component BECs [ 14] (see also Ref. [4]), by 
construction, this mode captures the out-of-phase motion of 
the DB-soliton pair Furthermore, by symmetry, the in-phase 
oscillation of the DB-soliton pair in the trap will be performed 
with the frequency Wosc- These two characteristic frequencies 
coincide with the eigenfrequencies of two anomalous modes 
of the BdG spectrum of the trapped DB-soliton pair (see, e.g., 
Refs. ||3,[3[3l for a relevant discussion of such modes). The 
rest of the spectrum will still be discrete (due to the presence 
of the harmonic trap — see also Sec. Ill), but in this case also 
collisions of these anomalous modes of the DB-solitons with 
the modes of the background may induce instabilities (see also 
Ref. H) —see below. 

We now turn to a systematic numerical investigation of the 
above features and of the multiple-DB-soliton states. At first. 
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FIG. 4: (Color online) The left and right columns correspond, respec- 
tively, to an in-phase and an out-of-phase dark-bright soliton pair in a 
harmonic trap with Q. — 0.1. The top row of panels depicts the pro- 
files of the DB-soliton pairs (solid blue lines and dashed green lines 
corresponding, respectively, to the dark and bright components) and 
the trapping potential (dashed-dotted red line). The middle row of 
panels depicts the spectral plane (Ar,Ai) rescaled by the trap fre- 
quency Q,. The bottom row of panels depicts the numerical (small 
stars in red) and the analytical (circles in blue) results for the equilib- 
rium distance between the solitons as a function of /12 ; the theoretical 
prediction is based on Eq. ( I33t . 



we consider the two-DB-soliton state in the trap, results for 
which are summarized in Figs. |4]and|5] both for the in-phase 
and the out-of-phase configurations. In particular, the top left 
and right panels of Fig. |4] show examples of an in-phase and 
an out-of-phase stationary DB-soliton pair, respectively (both 
for fii = 3/2 and /i2 = 1). The two middle panels illus- 
trate the corresponding spectral planes, showcasing the linear 
stability of these configurations. The bottom panels of the fig- 
ure show the equilibrium positions of the soliton centers. In 
the in-phase case (bottom left panel), it is observed that larger 
chemical potential (number of atoms) in the second compo- 
nent leads to stronger repulsion and, hence, larger distance 
from the trap center In the out-of-phase case (bottom light 
panel), we observe a similar effect but in the reverse direction 
(due to the attraction of the out-of-phase bright-soliton com- 
ponents) for smaller values of the chemical potential. Notice 
that in both cases a good agreement is observed between the 
numerically observed equilibrium separations and the theoret- 
ically predicted ones from Eq. ( |33] l. 

To study the validity of Eq. (|34] | — pertinent to small- 
amplitude oscillations around the fixed points — we show in 
Fig. |5] the eigenvalues A of the excitation spectrum [both for 
the in-phase (left column) and for the out-of-phase (right col- 
umn) cases] as functions of /i2. The imaginary and real part. 
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FIG. 5: (Color online) The left and right columns of panels corre- 
spond, respectively, to an in-phase and an out-of-phase dark-bright 
soliton pair in a harmonic trap with Q = 0.1. Shown are the imagi- 
nary (top row of panels) and the real (bottom row of panels) parts of 
the eigenvalues as functions of for fii — 3/2. In the top panels, 
the theoretical predictions for the eigenfrequencies of the anomalous 
modes of the system, pertaining to the in- and out-of-phase oscilla- 
tions of the DB-solitons are depicted by dashed (red) lines. Notice 
that collisions of modes (eigenvalue crossings) observed in the top 
panels indicate the emergence of instability windows observed in the 
bottom panels. The instabilities are of the Hamiltonian-Hopf type 
and result in the emergence of eigenvalue quartets. 



Xi and Xr, of the respective eigenvalues, normalized over the 
trap strength ft, are respectively shown in the top and bottom 
panels of Fig. |5] In the top panels, it is straightforward to 
compare the analytical result of Eq. (|34t with the BdG result, 
namely the second anomalous mode of the spectrum, corre- 
sponding to the out-of-phase oscillations of the DB-soliton 
pair. Once again, good agreement is observed between the 
two; the differences may be partially attributed to the "inter- 
action" (i.e., collisions) of these modes with other modes of 
the BdG spectrum. It is clear from the comparison of the cor- 
responding columns that there exist narrow instability win- 
dows, arising due to the crossing of the anomalous mode(s) of 
the DB-soliton pair with eigenmodes of the background of the 
two-component system. These instabilities arise in the form of 
Hamiltonian-Hopf bifurcations iIstIi through the emergence of 
quartets of complex eigenvalues resulting from the collision of 
two pairs. The growth rates of the pertinent oscillatory insta- 
bilities are fairly small (i.e., the instabilities are weak) in both 
the in- and out-of-phase cases; it should be noted, however, 
that in the latter case, the formation of the quartets appears to 
be occurring in very narrow intervals. 

Naturally, the above considerations can also be generalized 
to three- or more DB-solitons, although the analytical calcu- 
lations become increasingly more tedious; again, as we will 
show below, in-phase or out-of-phase configurations are pos- 
sible in the presence of the trap. Pertinent examples, showing 
two different three-DB-soliton configurations, are illustrated 
in Fig. |6] In particular, the first column in the figure corre- 
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FIG. 6: (Color online) The left and right columns of panels corre- 
spond, respectively, to an in-phase and an out-of-phase three-DB- 
soliton configurations. The top row of panels depicts the respective 
stationary states, for /ii — 3/2, /i2 = 1 and Q = 0.1; solid (blue) 
lines depict the dark-soliton components, dashed (green) lines the 
bright ones, while the dashed-dotted (red) line shows the harmonic 
trap. The second row of panels depicts the spectral planes for the 
above stationary states, while the third and fourth rows of panels are 
equivalent to those of Fig.|5] but for the three-DB-soliton configura- 
tions. 



sponds to the in-phase three-DB-soliton state, while the sec- 
ond column corresponds to the out-of-phase variant thereof. 
In the case under consideration, there exist narrow parametric 
intervals of dynamical instability, which are narrower for the 
out-of-phase case (as in the case of the two-DB-soliton states). 
We should mention, in passing, that the dynamics of two- 
and three-DB soliton configurations was recently studied in 
Ref . ll26ll ; our study complements the latter by yielding analyt- 
ical approximations and a numerical continuation/bifurcation 
approach towards such states. 
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VI. CONCLUSIONS AND DISCUSSION 

In the present work, we have studied multiple quasi- 
one-dimensional dark-bright (DB) solitons in atomic Bose- 
Einstein condensates. Our theoretical results were motivated 
and supported by the experimental evidence of the formation 
of DB-soliton clusters in a two-component, elongated rubid- 
ium condensate, confined in a harmonic trap. The theoret- 
ical analysis was based on the study of two coupled, one- 
dimensional Gross-Pitaevskii equations. 

Starting from the case of a homogeneous condensate (i.e., 
in the absence of a trapping potential), we have employed a 
Hamiltonian perturbation theory to analyze the interaction be- 
tween two DB-solitons. Assuming that the DB-solitons are 
of low velocity and sufficiently far from each other, we have 
found approximate expressions for the interaction forces be- 
tween the same or different soliton components. This way, we 
derived a classical equation of motion for the center of mass 
of the DB-soliton pair, and revealed the role of the phase- 
difference between the bright-soUton components: we have 
shown, in particular, that the repulsion between the dark soli- 
ton components may be counter-balanced by the attraction be- 
tween out-of-phase bright components, thus inducing the ex- 
istence of stationary DB-soliton pairs even in the case when 
the external trapping potential is absent. We have found the 
equilibrium distance between the two DB solitons that com- 
pose the stationary DB-soliton pair, with the semi-analytical 
result being in excellent agreement with the relevant numer- 
ical one. Additionally, we have demonstrated the linear sta- 
bility of these stationary DB-soliton pairs by means of ana- 
lytical and numerical techniques [the latter were based on a 
Bogoliubov-de Gennes analysis]. It was shown that the ana- 
lytical result for the oscillation frequency of small-amplitude 
perturbations around the equilibrium distance is in excellent 
agreement with the pertinent eigenvalue characterizing the ex- 
citation spectrum of the DB-soliton pair 

We have then studied multiple-DB-solitons in the trap. In 
this case, we have employed a simple physical picture, where 
the total force acting on the DB-solitons was decomposed to 
an interaction force (derived in the homogeneous case) and a 
restoring force induced by the trapping potential; the relevant 
characteristic frequency associated with the latter was the os- 
cillation frequency of a single-DB-soliton in the trap (which 
was found to coincide with the pertinent anomalous-mode 
eigenvalue of the single DB soliton system). Following this 
approach, we were able to find stationary in-trap DB-soliton 
pairs even in the case where the bright-soliton components 
were repelling each other: in this case, the trap-induced re stor- 



ing force was able to counter-balance the repulsive forces be- 
tween the dark- and the bright-soliton components. The semi- 
analytical results for the equilibrium distance and the oscil- 
lation frequencies (for the in- and out-of-phase cases) were 
again found to be in very good agreement with respective 
numerical results, including the anomalous modes' eigenfre- 
quencies pertaining to the in- and out-of-phase motion of soli- 
tons. The stability analysis of the DB-solitons in the trap 
indicated the possibility of the existence of unstable modes 
through Hamiltonian-Hopf instability quartets, although the 
latter would typically only arise over narrow parametric inter- 
vals — and with rather weak instability growth rates. Results 
pertaining to three-DB-solitons in the trap were presented as 
well; the main features of these states were found to be quali- 
tatively similar to the ones of the DB-soliton pairs. 

Coming back to our experimental findings, we should note 
the following. Given the nature of the available initial con- 
ditions, our experimental observations were not able to iden- 
tify, in a straightforward way, genuinely stationary DB-soliton 
complexes and/or to identify precisely their internal modes. 
Nevertheless, the frequent and persistent occurrence of DB- 
soliton clusters in the experiment is highly indicative of the 
robustness of such "DB-soliton molecules". This is in tune 
with our existence and stability results. 

It would be particularly interesting to further explore the 
dynamics of multiple-DB-soliton complexes, and potentially 
the formation of "DB-soliton gases" comprising such inter- 
acting atomic constituents. Deriving Toda-lattice-type equa- 
tions describing such gases, and identifying their stationary 
states, excitations and (mesoscopic}solitons (as in the case 
of single-component dark solitons Il35ll ). would be challenges 
for future work. Another possibility is to extend the present 
considerations to the vortex-bright solitons found in Ref. [28(]. 
There, it would be relevant to identify whether molecular 
states consisting of two- or of three-vortex-bright solitons can 
be constructed, and whether the relative phases of and tt be- 
tween the bright components can still yield different station- 
ary states. Relevant studies are presently in progress. 
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